Inference of Time-Evolving Coupled Dynamical Systems in the Presence of Noise 
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A new method is introduced for analysis of interactions between time-dependent coupled oscillators, based 
on the signals they generate. It distinguishes unsynchronized dynamics from noise-induced phase slips, and 
enables the evolution of the coupling functions and other parameters to be followed. It is based on phase 
dynamics, with Bayesian inference of the time-evolving parameters achieved by shaping the prior densities to 
incorporate knowledge of previous samples. The method is tested numerically and applied to reveal and quantify 
the time-varying nature of cardiorespiratory interactions. 
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The common assumption that a dynamical system under 
study is isolated and autonomous is never rigorously true. 
Furthermore, it is often a poor approximation, because the 
inevitable external influences may be too strong to ignore. 
For an oscillatory system, they can e.g. modify its natural 
frequency and/or amplitude. Much effort has therefore been 
made to understand non-autonomous oscillators driven from 
equilibrium by a variety of external forcings. A more difficult 
problem is faced where two or more interacting oscillatory 
systems are subject to external deterministic influences, a sce- 
nario that often arises in practice, e.g. in physiology includ- 
ing cellular dynamics, blood circulation, and brain dynamics. 
In such cases, the interacting systems (e.g. cardio-respiratory) 
are influenced by other oscillatory processes as well as by 
noise. Similarly, interactions at the intercellular level HI] and 
between subcellular components J2l are crucial to multicellu- 
lar organisms. Evaluation of the interactions by analysis of 
physiological signals and references therein) has proved 
useful in relation to a diversity of different diseases. 

Granger causality |0, Etl an d transfer entropy [0 0] have 
brought insight into the functional connectivity of systems, 
especially in neuroscience. Based on autoregressive and 
information-theoretic approaches to data-driven causal infer- 
ence, these methods focus on the statistical properties of the 
time series by measuring the extent to which the individual 
components exchange information. However, these methods 
are designed to infer effect, not mechanism. In contrast, we 
consider here complex interacting systems that are oscillatory 
and subject to noise, and extract their dynamical properties. 

Several questions immediately arise in relation to the dy- 
namics of coupled systems. Does the external influence alter 
their natural frequencies or amplitudes? Are they synchro- 
nized, or do they exhibit finite coherence? If synchronized, 
is it continuous or only for some of the time? Measurements 
may be relatively straightforward, using modern sensors and 
digital signal acquisition equipment, but how are the resultant 
signals to be analysed to reveal the characteristics of the origi- 
nating systems? To date, this inverse problem has no solution. 

Earlier work on coupled oscillators emphasized the detec- 
tion of synchronization lit-ITDl. and quantifying the couplings 
and directionality of influence between the oscillators 1112 - 



1511 . The inference of an underlying phase model enabled ex- 
traction of the phase-resetting curves, interactions and struc- 
tures of networks IU6H21I1 . However, these techniques inferred 
neither the noise dynamics nor the parameters characterizing 
the noise. In a quite separate line of development, Bayesian 
inference I22M27I1 has opened the door to the analysis of noisy 
time-evolving phase dynamics. 

In this Letter we introduce a new method that (a) encom- 
passes time-variable dynamics, (b) detects synchronization 
where it exists, and (c) determines the inter-oscillator coupling 
functions regardless of whether or not they are time varying. 
By reconstructing the dynamics in terms of a set of base func- 
tions, we evaluate the probability that they are driven by a set 
of equations that are intrinsically synchronized, distinguish- 
ing phase-slips of dynamical origin from those attributable 
to noise. The Bayesian probability lying at the core of the 
method is itself time-dependent via the prior probability as a 
time-dependent informational process. Thus relatively small 
windows can provide good time-resolved inference. 

When two noisy, weakly-interacting, ./V-dimensional, self- 
sustained oscillators synchronize |28|] . their motion is de- 
scribed by their phase dynamics: 



4>i =uj t + fi (&) + gi(<f>i, fa) + £i(t) 



(1) 



leaving other coordinates expressed as functions of the phase: 
rj = rj(0j). £ is a two-dimensional noise, usually assumed 
Gaussian and white, (£i(t)£j(T)) = S(t — r)Eij and which 
may, or may not be, spatially correlated. Noise can induce 
phase slips in a system that would be synchronized in the 
noise-free limit, so evaluation of synchronization needs pre- 
cise inference of fi and gt, and of the noise matrix Ey. The 
systems' periodic nature suggests periodic base-functions, 
whence the use of Fourier terms for the decomposition: 
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Assuming that the dynamics is adequately described by a fi- 
nite number K of Fourier terms, we can rewrite the phase 



2 



dynamics of (Q]) as a finite sum of base functions: 



K 



J2 cj? , S Ilfc (0i J a ) + £i(t). 



(3) 
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where I = 1,2, $i, = $2. 
and c^, are the if most important Fourier components. 

In order to reconstruct the parameters of (H) we exploit the 
approach already presented in J25 26 1 assuming that a 2- 
dimensional time-series of observational data X = {<fii n = 
<M^n)} = wft) is provided, and that the unknown model 
parameters M. = {c? , i^y} are to be inferred. 

In Bayesian statistics a given prior density p pi i or (A4) 
that encloses expert knowledge of the unknown param- 
eters (based on previous observations) and the likeli- 
hood function £(X\A4), the probability density to observe 
{<t>i,n{t)} given choice M. of the dynamical model, are 
used to calculate the so-called posterior density px{M.\X) 
of the unknown parameters M., conditioned on observa- 
tions, by application of Bayes' theorem px{M.\X) — 
£(X\M) Ppnoi (M)/ J £(X\M) Ppnm -(M)dM. 

For independent white Gaussian noise sources, and in 
the mid-point approximation where = ^ l ' n+1 h and 
't'ln = {'t'l-n + 0;,n+i)/2, the likelihood is given by a prod- 



uct over n of the probability of observing 1 
The negative log-likelihood function S = 



_„_!_! at each time. 

'hi£(X\M) is 
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S=-ln\E\ + - ^(c fc + 
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with implicit summation over repeated indices k,l,i,j. The 
log-likelihood is a function of the Fourier coefficients of the 
phases. Hence for a multivariate prior probability, the pos- 
terior probability is a multivariate normal distribution. From 



I 25L 12611 . and assuming such a distribution as a prior for pa- 
rameters Cj. , with mean c, and covariances H _1 pi i or , the sta- 
tionary point of S is calculated recursively from: 



/t d$t,fc(0.,n) 

2 d<j>i ' 

= Sprier^ + ^<M<J ^>(<£*J, 



(4) 



with implicit summation over 11 = 1, . . . , N and over repeated 
indices k,l,i,j,w. The mean parameter vector of the posterior 



isthenc« = (S-inr 
"flat" prior can be used as the initial limit of an infinitely large 
normal distribution, by setting H pr i or = and c pr j or = 0. The 
multivariate probability Afx(c\, c, S) for the given time series 
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We note that a noninformative 



X explicitly defines the probability density of each parameter 
set of the dynamical system. 

When the sequential data come from a stream of measure- 
ments providing multiple blocks of information, one applies 
(0]i to each block. If the system is known to be non-time- 
varying, then the posterior density of each block is taken as 
the prior of the next one. Thus, the uncertainties in the param- 
eters steadily decrease with time as more data are included. 

If the system has time dependence, however, the method 
of propagating knowledge about the state of parameters obvi- 
ously has to be refined. Our framework prescribes the prior 
to be multinormal, so we synthesize our knowledge into a 
squared symmetric positive definite matrix. We assume that 
the probability of each parameter diffuses normally with a 
known diffusion matrix Sdiff. Thus, the probability density of 
the parameters is the convolution of two normal multivariate 



distributions, E post and S diff : 



^post 



The particular form of Ediff describes which part of the dy- 
namical fields defining the oscillators has changed, and the 
size of the change. In general (Sditf)i,j = pij<Ji(Tj, where o~i 
is the standard deviation (SD) of the diffusion of c; in the time 
window t w , and pij is the correlation between the change in 
the parameters c, and Cj . We will consider a particular exam- 
ple of S d iff: we assume there is no change of correlation be- 
tween parameters = <Sy) and that each SD cr, is a known 
fraction of the relevant parameter, cr, = p w Ci, where p w indi- 
cates that the parameter p refers to a window of length t w . 

The probability of synchronized dynamics is estimated by 
sampling the posterior and evaluating its overlap with the 
Arnold tongue border: p sync = / s(c) Mx(c\c, 5) dc, where 
s(c) = {1, 0} defines whether the parameter set c is inside or 
outside the synchronization region. For motion on the torus 
T 2 defined by the toroidal coordinate C(0i(i), 4>2(t)), and the 
polar coordinate ip(t), we consider a Poincare section defined 
by C = and assume that d^(t) / dt\^—o > for any ip. Thus 
the direction of motion along the toroidal coordinate is the 
same for every point of the section, which we would like to 
follow in order to check whether there is a periodic orbit. If 
so, and if its winding number is zero, then the system is syn- 
chronized; and there must be at least one other periodic orbit 
with one of them being stable and the other unstable. 

Solution of the dynamical system over the torus yields a 
map M: [0, 2ir] — > [0, 2ir] that defines, for each ip n on the 
Poincare section, the next phase %p n +i after one period of the 
toroidal coordinate: ip n +i = M(ip n ). The map M is contin- 
uous, periodic, and has two fixed points (one stable and one 
unstable) if and only if there are two periodic orbits for the 
dynamical system, i.e. synchronization is verified if ip e exists 
such that ip e = M(ip e ) and \dM(ip e )/dijj\ < 1. To calcu- 
late s(c) for any of the sampled parameter sets, we: (i) fix an 
arbitrary ( and, for any tpi, integrate (0 numerically for one 
cycle of the toroidal coordinate, obtaining the mapped point 
M(ipi); (ii) by finite difference evaluation of dM/dtp employ 
a modified version of Newton's root-finding method to find 
the occurrence (if any) of ip such that M (ip) = ip. If there is 
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a root, s(c) = 1 is returned, otherwise s(c) = is returned. 

From the inferred parameters of the base functions fi(4>i), 
gi((f>i, 4>j), we can reconstruct the specific functional form of 
the coupling functions qi(<f>i,<f>j). The novel advantage of 
this framework is that it allows reconstruction of the time- 
variability and evolution of such coupling functions. Simple 
normalization of the inferred coupling parameters yields the 
inter-oscillator coupling strengths, and thence the directional- 
ity index Il2l - ll5ll . If D € (0, 1] the first oscillator drives the 
second (1 — > 2), or if D e [—1,0) otherwise. Note that, al- 
though our discussion relates to two oscillators, Eqs. ([I}-© 
are also applicable to a network of oscillators. For expanded 
discussion, technical description and software codes see ll37ll . 

As a demonstration of how the synchronization detection 
works, we simulated numerically a pair of coupled, noisy, 
phase oscillators (Q~|i. Bayesian inference followed by exam- 
ination of the constructed map M(tp e ) showed that our ap- 
proach successfully distinguishes synchronized (s(c) = 1) 
from unsynchronized dynamics (s(c) = 0), i.e. whether the 
root M(ip e ) = ipe exists or not. To demonstrate the nov- 
elty of our method we consider the characteristic case illus- 
trated in Fig.Q] The parameters were such that the oscillators 
were only just inside the Arnold tongue so that, for moderate 
noise, phase slips occurred, as shown schematically in Fig. 
Ola). The application of earlier methods based on the statis- 
tics of the phase difference lli UIoll suggests that the oscillators 
are not synchronized. In contrast, our new technique shows 
that the oscillators are intrinsically synchronized as shown in 
Fig. 0I C ) : the phase slips are attributable purely to noise (the 
intensity of which is inferred in matrix Ei j), and not to deter- 
ministic interactions between the oscillators. 

To see how the new method can also follow time-variations 
of the parameters, coupling functions and synchronization, we 
take as an example two coupled noisy Poincare oscillators: 

X i = ~(\J X i + Vi ~ l ) X i ~ <^(%« + e i{ t ){ x j ~ X i) + &(*) 

tfi = -(yo? +vf - l )Vi + Ui(t)xi + ei(t)(yj - yi) + £i(t) 
with* = 1,2; j = 1,2; i^j. (5) 
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Figure 1: Synchronization discrimination for the coupled phase os- 
cillators {Q with: wi = 1.2, u 2 = 0.8, ei = 0.1, £2 = 0.35, 
fi{4>i) = 0, gi = sin(<^>2 - 4>i), 92 = sin($i - and noise 
strengths En = E22 = 2. (a) Schematic Arnold tongue to il- 
lustrate synchronization f29jl. (b) Phase difference, exhibiting two 
phase slips, (c) Map of M(tp e ) for (b) demonstrating that a root of 
M{ip) — ip exists, i.e. that the state is in fact synchronized. 












Figure 2: Extraction of time-varying parameters, synchronization 
and coupling functions from numerical data created by The 
plots show the results inferred for the numerical values of constants 
listed in the text. The frequency fi(t) and coupling £2(t) are in- 
dependently varied: (a) Wi(t) = Wi + Ai sin(wit); (b) E2(t) = 
£2 + A2 sm(u)2t). The dotted and full lines plot the parameters when 
the two oscillators are synchronized for part of the time (ei = 0.3), 
and not synchronized at all (ej = 0.1), respectively. The regions of 
synchronization, found by calculation of the synchronization index, 
are indicated by the gray shaded regions, (c)-(f) show the coupling 
functions qi(<j>i,(j>2) and 92(^1, 4>2) for time windows centered at 
different times: (c) and (d) at t = 350s; (e) and (f) at t = 1000s. The 
window length t w — 50s, and £1 =0.1 in both cases. Note the simi- 
larity in forms of (c) and (e), and of (d) and (f). The other parameters 
were: e 2 = 0.1, Ui = 2tt1, cj 2 = 2tt1.14, Ai = 0.2, A 2 = 0.13, 
uii = 2tt0.002, <Z>2 = 2tt0.0014 and noise E n = E22 = 0.1. The 
phases were estimated as <j>i = arctan(j/i/a;i). 



We consider bidirectional coupling (1-0-2), where the natu- 
ral frequency of the first oscillator, and its coupling strength 
to the second one, vary periodically. For e\ = 0.1 there is 
no synchronization: the time-varying parameters (fi(t) and 
e^it)) are accurately traced (full red lines of Fig. |2ja) and 
(b)). For a coupling of £\ = 0.3 the two oscillators will 
be synchronized for part of the time, resulting in intermit- 
tent synchronization. The time-variability of the parameters in 
the non-synchronized intervals is again determined correctly 
whereas, within the synchronized intervals, the inferred pa- 
rameters (dashed lines in (a), (b)) diverge from their true val- 
ues (full red curves). Within these synchronized intervals, all 
of the base functions are highly correlated, with values lying 
within the Arnold tongue. The latter was detected as the range 
for which s(c) = 1, grey-shaded in Figs. [2 a) and (b). 

The reconstructed sine-like functions 91(^1,^2) and 
52(01,02) are shown in Figs. |2jc) and (d) for the first and 
second oscillators, respectively. They describe the functional 
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Figure 3: Synchronization, directionality and coupling functions in 
the cardiorespiratory interaction, (a) Standard l:N synchrogram. 
(b) Synchronization index for ratios 1:4, 1:5 and 1:6, as indicated. 
The light-gray dotted line represents the mean, and the dark-gray 
dashed line the mean +2 SD, of synchronization indices calculated 
from 100 surrogate [f30h realizations, (c) The time-varying respira- 
tion frequency (note the downward ramp due to pacing). The gray 
areas in (c) represent ±2 SD from the mean value, (d) Directional- 
ity index (full curve); the light-gray dotted line represents the mean 
directionality index calculated from 100 surrogate realizations, and 
the dark-gray dashed line represents the mean +2SD. (e)-(g) cou- 
pling functions qi (<f>i, cj>2) calculated at different times, as indicated 
by the grey arrows. 



form of the interactions between the two Poincare systems ©. 
The results suggest that the form of the coupling functions 
does not evolve with time: qi and q2, evaluated for later time 
segments, are presented in Figs.|2](e) and (f) respectively. By 
comparison of Figs. |2jc) and (e), or of Figs. [2jd) and (f), we 
see that the coupling functions did not change qualitatively, 
even though there were time-varying parameters and weak ef- 
fects from the noise. 

It is well known that modulation and time-variations tend 
to affect synchronization between biological oscillators 
3211 . Hence the need for a technique able, not only to identify 
time-varying dynamics, but also to evaluate measures of inter- 
action, e.g. synchronization, directionality and coupling func- 
tions. To demonstrate the method on real biological data, we 
analyzed cardio-respiratory measurements from resting hu- 
man subjects whose paced respiration was ramped down with 
decreasing frequency. The instantaneous cardiac phase was 
estimated by wavelet synchrosqueezed decomposition i33ll of 
the ECG signal. Similarly, the respiratory protophase was ex- 



tracted from the CO2 concentration signal, followed by trans- 
formation lfl9ll to the phase. The results are shown in Fig. 
m First, just for comparison, the corresponding synchrogram 
12811 of the same data is presented in (a). The time-variation 
of the respiration frequency is clearly evident in (c). By nor- 
malizing the inferred coupling parameters, we determined the 
net directionality of the interactions. Fig. [3}d) suggests that 
the degree of directionality is time- varying; the analyses con- 
firm that respiration-to-heart is dominant I13l ll2l - ll5ll . even for 
non-paced respiration (not shown). The set of inferred pa- 
rameters and how they are correlated can be used to deter- 
mine whether cardiorespiratory synchronization exists and, if 
so, in what ratio. Fig. [3jb) shows transitions from the non- 
synchronized to the synchronized state, in ratios 1:4 to 1:5 
to 1:6, as the ramp progressed. The cardio-respiratory cou- 
pling function, evaluated for three different time windows is 
presented in Fig. [3e)-(g). Note that the interactions are now 
described by complex functions whose form changes qualita- 
tively over time - cf. Fig. |3]e) with (f) and (g). This implies 
that, in contrast to many systems with time-invariant coupling 
functions (e.g. Fig.[2jc)-(f) or |34j-|36J]), the functional rela- 
tions for the interactions of an open (biological) system can 
themselves be time-varying processes. By analyzing consec- 
utive time windows, we can even follow the time-evolution of 
the coupling functions - cf. the similarities i.e. evolution of 
Figs. [3 f) and (g). It is important to note that the variability in 
form of the coupling function can cause synchronization tran- 
sitions. This variability is not caused by the time-varying res- 
piration frequency (which is decomposed separately). We also 
observed time-evolution of the coupling functions for sponta- 
neous (non-paced) breathing. 

In summary, our new method for inference of phase dy- 
namics enables the evolution of a system to be tracked con- 
tinuously. Unlike earlier methods that only detect the occur- 
rence of transitions to/from synchronization, it reveals details 
of the phase dynamics, describing the inherent nature of the 
transitions and simultaneously deducing the characteristics of 
the noise that stimulated them. We have identified the time- 
varying nature of the functions that characterize interactions 
between open oscillatory systems. The cardio-respiratory 
analysis demonstrated that not only the parameters, but also 
the functional relationships, can be time-varying, and the new 
technique follows their evolution effectively. This novel fa- 
cility immediately invites many new questions, e.g. the func- 
tional forms between which the couplings vary, their frequen- 
cies of variation, how their variation affects synchronization 
transitions, and whether there is periodicity or a causal rela- 
tionship waiting to be identified and understood. Thus a whole 
new area of investigation has become accessible. 

Our grateful thanks are due to Dwain Eckberg for providing 
the data for Fig. [3] and to M. Arrayas, M.I. Dykman, M. Hor- 
vat, D. Iatsenko, RE. Kloeden, D.G. Luchinsky, R. Mannella, 
S. Petkoski, and M.G. Rosenblum, for valuable discussions. 
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